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Abstract 

Fourier acceleration has been successfully applied to the simulation of lattice field theories for 
more than a decade. In this paper, wc extend the method to the dynamics of discrete particles 
moving in continuum. Although our method is based on a mapping of the particles' dynamics to 
a regular grid so that discrete Fourier transforms may be taken, it should be emphasized that the 
introduction of the grid is a purely algorithmic device and that no smoothing, coarse-graining or 
mean-field approximations are made. The method thus can be applied to the equations of motion 
of molecular dynamics (MD), or its Langevin or Brownian variants. For example, in Langevin 
MD simulations our acceleration technique permits a straightforward spectral decomposition 
of forces so that the long-wavelength modes are integrated with a longer time step, thereby 
reducing the time required to reach equilibrium or to decorrelate the system in equilibrium. 
Speedup factors of up to 30 are observed relative to pure (unaccelerated) Langevin MD. As with 
acceleration of critical lattice models, even further gains relative to the unaccelerated method 
are expected for larger systems. Preliminary results for Fourier-accelerated molecular dynamics 
are presented in order to illustrate the basic concepts. Possible extensions of the method and 
further lines of research are discussed. 
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1 Introduction 



Molecular dynamics (MD) simulations play an important role in our fundamental understanding of 
the kinetics of molecular systems and provide a powerful tool for modeling a wide variety of materials 
including biomolecules. Although MD simulations have benefited tremendously from advances in 
high-performance computing, they suffer from the limitation arising from the numerical stiffness 
inherent in Newton's equations. The result is that MD studies are generally restricted to short 
intervals of real time, from nanoseconds up to a few microseconds, even with heroic computational 
efforts. To overcome this difficulty, there is a growing effort to develop accelerated MD algorithms. 
(See, for example, ) 

In contrast to molecular (or other discrete particle) systems with a Lagrangian data repre- 
sentation, there is a considerable variety of acceleration algorithms available for continuum field 
theories approximated on a regular grid or lattice. For example, grid-based simulations have made 
substantial progress with the advent of cluster Monte Carlo methods ||^, ^, Fourier acceleration Q, 
and multi-grid iterative solvers |1C]. Because the bulk properties of large aggregates of molecules 
can often be described by continuum mechanics, it is intuitively appealing that a corresponding 
method should apply in the molecular (or particulate) framework. Indeed, making this connection 
between the molecular and continuum scales is a central goal for multi-scale modeling projects. In 
this paper we show how one such continuum tool, namely Fourier acceleration, can be applied to 
Langevin MD without introducing any coarse-graining or mean- field approximations. The basic 
ingredient of the method is an exact mapping of the original particulate system onto a regular 
lattice of displacement fields. Although this mapping may prove useful in a broader context, we 
restrict our attention to Fourier acceleration of the Langevin equations for MD. 

The idea of introducing a regular grid into MD is not new; grid-based recursive multipole 
expansions [11 1, for example, have been used for more than a decade to rapidly compute Coulomb 
interactions. More recently, hybrid atomistic-continuum techniques, such as the quasi-continuum 
method use finite-element techniques to bridge microscopic and macroscopic length scales. 
Most applications of spectral methods to molecular systems, however, have been confined to the 
analysis of data (for example, structure and response functions). 

By contrast, our procedure uses spectral analysis to modify and accelerate the dynamical evo- 
lution of the molecular system. Unique to this approach is the mapping of the actual position 
coordinates to a grid, and the ability to invert the mapping to displace the original off-lattice 
molecular coordinates. The introduction of the grid is a purely algorithmic device and is not tan- 
tamount to a coarse-graining or mean- field approximation of any kind; that is, the accelerated 
dynamics are still those of discrete particles. The result is a new, accelerated, stochastic dynamics 
that is significantly faster than standard Langevin MD, but still exactly preserves the equilibrium 
distribution. The fundamental tradeoff associated with this approach is that of speed versus faith- 
fulness to the essential kinetics. Both of these desiderata are clearly specific to the system being 
studied and the phenomena that the model should faithfully represent. 

The organization of this paper is as follows. Section |^ describes our procedure for mapping 
the particulate system to a regular lattice. This mapping is a prerequisite to the application of a 
Fourier-mode decomposition. In Section |3| we outline the Fourier-accelerated Langevin dynamics on 
the grid. We demonstrate the method in Section ^ by applying it to a cj)'^ model at its critical point. 
We then describe how to apply Fourier acceleration (FA) in conjunction with the lattice mapping in 
Section ^. As an example, in Section ^, we apply the method to the Langevin dynamics of a simple 
Lennard-Jones fluid. Extensions of the Fourier-accelerated molecular dynamics (FAMD) method 
and additional applications are discussed in Section ^. 
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2 Particle-to-Grid Mapping Procedure 



There are many ways by which a molecular system can be transformed from (off-lattice) particle 
coordinates to a fiducial grid. Each has its advantages and disadvantages, depending upon the aim 
of the transformation. In this section we discuss one method that has proven to be particularly 
useful. 

In one dimension, the simplest mapping procedure procedure is to sort the particles by their 
position coordinate. Each particle i is given a permuted label n{i) so that n(z) < n{j) if Xi < Xj. 
Whereas the i and j indices are arbitrary labels, devoid of physical significance, the permuted labels 
are based on the sequence of particle positions and hence may be thought of as lying on a grid 
with some physical meaning. Because n{i) is a permutation, it has inverse function i{n) such that 
n(z(-)) = -. We now transform to new coordinates by the prescription Xn = Xi[n)- This mapping 
makes it possible to directly Fourier transform the new position coordinates, 

Xk = \^Y.^ne'^''- (2-1) 

n 

Note that this new spatial representation contains precisely the same amount of information as the 
original data. Also note that because this mapping is merely a geometrically motivated relabeling, 
all attributes, such as mass m,, are automatically transfered to the new representation. 

The multidimensional generalization of this method is not so straightforward. The problem of 
sorting the particles in more than one dimension is not well defined. There is, however, a very good 
and efficient approximation used by numerical analysts for load-balancing graphs on multiprocessor 
architectures, known as Recursive Coordinate Bisection (RGB) |T^. To see how this algorithm 
works, consider a two-dimensional square domain containing N = L"^ particles, where L is a power 
of 2. We first introduce a two-index label i = {ix,iy) for the particles, where 

ix{i) = i mod L (2.2) 
iyii) = {i-ix)/L, (2.3) 

so that 



z(i) = ix + {iy-l)L. (2.4) 

Thus the transformation from one-index labels i to two-index labels i is a bijection. We can label 
the particles' coordinates as Xj = {xi,yi), or equivalently as Xi = {xi,yi) = (xj(i), yj(i)). as with the 
i labels for the one-dimensional case, these labels (both i and i) are assigned arbitrarily and devoid 
of physical content. 

Whereas it is difficult to see how to order the one-index labels i, the RGB method provides 
a straightforward prescription for permuting the two-index labels i into a new set of two-index 
labels n(i) that are based on the particles' positions. Again, this function is a permutation, so it 
has an inverse i(n), such that n(i(-)) = •. The n's may reasonably be taken to lie on a regular 
two-dimensional grid, and hence provide a set of independent variables with respect to which the 
new coordinates Xn = Xj^^) = Xj^j^n)) can be Fourier transformed, just as in the one-dimensional 
example. 

To accomplish this, the physical domain is first divided into left- and right-hand portions with 
equal numbers {N/2) of particles, by sorting the particles on their x coordinates. The half with 
smaller x coordinates will have 1 < < L/2 and the half with larger x coordinates will have 
L/2 < Ux < L. In binary notation, this labeling sets the most significant bit of the Ux index to 
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Figure 1: Grid Mapping for a two-dimensional Lennard Jones fluid by Recursive Coordinate Bisec- 
tion. 

zero/one for the left /right halves. Next, we sort each set of N/2 particles by their y coordinates to 
obtain 4 sets of N/i particles, likewise setting the most significant bit in the Uy. This procedure 
is then applied recursively to each of the 4 boxes with iV/4 particles, maintaining the alternation 
between the x and y axes. For systems of relatively uniform density, the resulting fiducial grid leads 
to a remarkably regular and local particle-labelling scheme. RGB is an order A^logA^ algorithm 
with many obvious similarities to fast Fourier transforms. 



3 Fourier Accelerated Langevin Dynamics 

To demonstrate how Fourier Acceleration works, we consider in detail a simple (discrete-time) 
Langevin dynamics. The Langevin equation of motion for a system of N particles is 

Xi{t + At) = Xi(t) + 1^ {Atf + Mt)At, (3.1) 
zrrii 

where the A'" momenta are Gaussian random variables (pi(t)pj(t')) = ^kBTrriiSijSt^t'^- It is well 
known that this dynamics (in the limit of vanishing time step) samples the canonical-ensemble 
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Boltzmann-Gibbs equilibrium distribution function, 



P(xi,pi) = — exp 



Irrii 



(3.2) 



where /3 = l/ksT and the force is fi = —dV/dxi. 

For the moment, we set this result aside and consider lattice field problems for which Batrouni 
et al. have shown how to accelerate the approach to equilibrium in Fourier space. For example, 
consider fields (/>x(i) on a uniform grid with sites x, obeying the equilibrium distribution. 



P( 



1 



,{t + At) = Mt) - K^§^ {Atf + 7i^{t)At , 



(3.3) 



with action S. This distribution is a fixed point of the discrete Langevin dynamics (as At — > 0), 

(3.4) 



where T/x(i) are Gaussian random fields. This Markov process, however, is not the only one which 
drives the system to the equilibrium in Eq. (|3.3| ). Indeed, the local dynamics of Eq. ( |3.4| ) generally 
exhibits long autocorrelation times near critical points. Batrouni et al. |^] have shown how to 
accelerate such grid-based Langevin equations using a Fourier decomposition of the dynamics. The 
Fourier Acceleration (FA) method depends on the simple observation that any mobility (or inverse 
mass) matrix may be introduced by the substitutions, K — > -ftTx.x' and (??x(i)'?x(^)) — K^^^i5{t' ,t)\, 
without upsetting the equilibrium distribution of the fields. One such choice is a matrix K which 
is diagonal in Fourier space. This choice leads to acceleration if the time steps for the slow (low 
wavenumber) modes are amplified. 



,{t + At) = ^^{t) + 



-k{\i)T 



dS 



{Aty + Vi^(k) ??k At 



(3.5) 



where J- represents a Fourier transform. A simple substitution of the field with the position 
of a particle Xi would allow us to use Fourier methods. The mappings of Sec. |2| provide that 
substitution. 



4 04 Model 

To illustrate the above procedure we apply the FA technique to the model at the critical point 
in two dimensions. It has been shown by Batrouni et al. that critical slowing down is completely 
eliminated by FA in a purely Gaussian model. Of course, in that case, the modes completely 
decouple in momentum space and each can be integrated independently. For a nonlinear model 
with mode-coupling, it is not guaranteed that FA will work at all. It is also not clear a priori 
what the optimal choice of the mass matrix should be that will most rapidly drive the system to 
equilibrium or decorrelate the system once in equilibrium. 

To gain experience in selecting the mass matrix for FAMD, we first studied a simpler system, 
namely the (j)'^ model in two dimensions at criticality. This model provides a qualitative (and in 
some cases quantitative) description of a displacive phase transition. The investigation of such 
transitions is one of our long-term goals. Surprisingly, the FA method applied to this system at its 
critical point has not been analyzed. However, Batrouni and Svetitsky have studied its application 
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to first-order phase transitions in a (j)'^ model and found a significant speedup of tunneling between 
minima |14]. 



The Hamiltonian is given by, 

mm) -- 



N 

E 

i=l 

N 

E 

i=l 



^'A. + 4<^. + 2E 

2d 



^ a2 I X ,4 

+ -^(Pi 



where 



2d-e . 



(4.1) 



(4.2) 



This system exhibits critical behavior along a line of critical parameters, x ^'^^ ^- We simulated 
the system at the critical point, x = 1-0, = 1.265, which was numerically determined previously 



by Toral and Chakrabarti [18|. 

We updated this system using the Fourier accelerated Langevin equation described above, 
namely, 

Ut + At) = Ut) + K{\<i)T (^-^) + ^UbT K(k) r?(k) , (4.3) 

where 77(k) represents the Fourier transformed Gaussian noise with {rf) = and {rf) = 1, and K(k) 
represents the mass matrix which gives us the desired acceleration. We chose the mass matrix to 
be the lattice propagator of the free theory. 



i^(k) 



Ad + m? 



,2 kfj_ 
2 



{Atf , 



(4.4) 



where the parameter m is expected to be order or 1/L for finite-size scaling Q. The value of 
this parameter was adjusted during trial simulations by setting m = c/L for different values of the 
constant c. We report the results for c = 4^/2. As a check, we repeated the simulations without 
Fourier acceleration using the pure Langevin update. 



4>i{t + At) = Ut) + —/('A) + ?7, 
where is a zero- mean unit-variance Gaussian random variable and the force term, /((/>), is. 



(4.5) 



5rL 



N 



2d 



xr = -E^>^-x^?-E 



(4.6) 



Finite-size scaling simulations were conducted for L x L system sizes where L = 2, 4, 8, 16, 
32, 64 using the FA Langevin update and for L = 2, 4 , 8, 16 for the pure Langevin case. A time 
step of At = 0.05 was used. Note that our time step corresponds to in Ref. Q and thus should 
give very little discretization error. We ran each system for time that is approximately 1000 times 
longer the correlation time estimated from trial simulations. The results are shown in Fig. ^. The 
normalized correlation functions, C{t) = {E{0)E{t)) / {E{0)E{0)) , were computed from the time 
series of the energy density. Correlation times r for each L were computed by fitting the region 
0.3 < C{t) < 0.6 to exp(— t/r). Error bars were estimated from the standard deviation of the 
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Value 


Error 


Value 


Error 


Value 


Error 


L = 


2 


0. 481 


± 0. 0159 


0. 496 


± 0. 0612 


0. 473 


± 0. 


0432 


L = 


4 


3. 174 


± 0. 0411 


3. 180 


± 0. 4019 


3. 291 


± 0. 


2998 


L = 


8 


14. 54 


± 0. 1576 


14. 65 


± 0. 9797 


14. 87 


± 0. 


5239 


L = 


16 


61. 16 


± 0. 4790 






62. 58 


± 1. 


6278 


L = 


32 


251. 5 


± 1. 6066 






260. 6 


± 3. 


5662 



Table 1: Mean energies and errors for each system size and update method. Errors are standard 
deviations from 10 blocked averages. 

values of r measured from five independent time series per system. Average energies and standard 
deviations of 10 blocked averages were measured for each system size and update algorithm. These 
results are listed in Table |l|. 

In Fig. ^ we compare the autocorrelation times for the pure Langevin update with those for the 
Fourier-Langevin case. There is clearly an acceleration. Whether the dynamical exponent z, which 
describes the growth of autocorrelation times by the scaling relation, r = L^, is actually different 
for Langevin and Fourier Acceleration is an interesting and open question, and would require more 
extensive computation than we have done to date. For practical simulations of systems far from 
criticality, the value of z is often not as important as the overall amplitude of the autocorrelation 
time. 

5 Fourier Accelerated Molecular Dynamics 

To apply these techniques to discrete particles with a Lagrangian data representation, we must 
introduce a Fourier transform of the position coordinates Xj(i). Clearly we cannot simply transform 
the Xj with respect to the particle labels i = 1^2, . . . N . As mentioned in Sec. ^ these labels are 
generally devoid of physical meaning. They have no natural relationship to the properties of the 
particles or to their spatial and/or temporal configuration. Hence, the first step is to map the 
particles onto a uniform spatial grid. 

The mapping scheme discussed in Sec. § suggests how to define appropriate grid coordinates. 
In the Index Method, the mapping, Xj Xn, is simply a relabling of the coordinates. (To 
be more explicit this notation for a 2D system is expanded into components: Xj = {xi,yj) and 

= {Xn^^ny,yna:,ny), where n = (rixjUy) is a two component integer vector. ) Consequently the 
Langevin dynamics is unaffected, 

Xn{t + At) = X„(t) + -^At^ + ^/k^ Vn(.t)At, (5.1) 

where we have introduced the normalized independent Gaussian noise with variance {rjrj) = 1. 
Because we have established a two-dimensional grid, we may now try to accelerate the dynamics 
simply by going to Fourier space as described above for a generic lattice field theory. The fields are 
now the position vectors of each particle. As we will demonstrate numerically in Sec. ^, this grid 
is indeed useful for a two-dimensional Lennard-Jones fluid. 
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Figure 2: Finite size scaling of the correlation time r with the linear dimension L for Langevin and 
FA Langevin simulations of 0^ theory; r is in units of At. 
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Wavenumber 


FAMD 


Lan 


^evin 




Value Error 


Value 


Error 


n = 1 


12000 ± 3000 


80000 


± 12000 


n = 2 


3300 ± 500 


22000 


± 4000 


n = 4 


1100 ± 130 


8600 


± 600 



Table 2: Correlation times in (integration time steps) for different wavelengths for = 16 particles. 

6 2D Lennard-Jones Fluid 

Motivated by the successful application of Fourier Acceleration to decorrelate lattice 0^ systems, 
we have tested its ability to reduce the autocorrelation time of a system of Lennard-Jones atoms 
in two dimensions using the Index Method. 

The Lennard-Jones interaction potential is given by 

^^^W=46(^^-^). (6.1) 

In our simulations we chose e = a = 1, a potential cutoff at 2.5(T, and worked at the liquid-vapor 
critical point, with temperature and density parameters T = = 0.47, pc = 0.35, respectively. 
Both pure Langevin MD and FAMD were tested for A = 16 and 64 particles with periodic boundary 
conditions. Each system was evolved on the order of 10^ integration steps with At = 0.005. 
This time step allowed us to accurately determine the critical thermodynamic quantities. The 
acceleration kernel we used in FAMD was identical in form to the one we applied to the (f)^ model, 
namely 

e(k) = ^ , (At)2 , (6.2) 

4E2.isin2(%)+^ 

where A^ is now the number of particles in the system. This should be compared to Eq. ( [4.4] ). 

We allowed the system to evolve for 10^ steps before statistics were taken. To compare the 
effectiveness of the Fourier acceleration, we examined the autocorrelation of various long- wavelength 
modes of the system. In particular, we looked at the circularly averaged time-autocorrelation of 
the cosine-transformed density. In D spatial dimensions, we write dk = k^~^dk dCl, where /c = |k| 
and dCt is the direction differential, so the cosine-transformed density is 

N 

p(k, t) = ^ cos(k • Xj). (6.3) 

i 

The autocorrelation that we measure is then 

A(^,r).I^^^l^%i±lMM^ (6.4) 

where k = 27rn/L, and L is the linear dimension of the system. 

In Tables | and |, we show the autocorrelation times for various modes in both Langevin and 
FAMD simulations. As seen from the tables, the FAMD dynamics is clearly more efficient at 
decorrelating long wavelength modes. A precise measure of the gain over standard Langevin MD 
was not possible, because standard Langevin MD has a very long correlation time. We therefore do 
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Wavenumber 



FAMD 



Langevin 



Value 



Error 



Value 



Error 



n = 1 
n = 2 
n = 4 



40000 ± 13000 
15000 ± 2000 
4200 ± 700 



1200000 ± 130000 
310000 ± 60000 
40000 ± 4000 



Table 3: Correlation times in (integration time steps) for different wavelengths for = 64 particles. 

not know exactly how much faster FAMD is. Moreover, whether or not there is simply a decrease 
in the amplitude of decorrelation time or an actual decrease in its algebraic form is not known. 
As with the precise determination of z for the (/)'^-model simulation, that will require considerably 
more computational effort which we postpone to future work. 

Finally, we limited our investigation to a maximum of only 64 particles to allow us to equilibrate 
the system at its critical point using Langevin dynamics. We expect, that gains over standard 
Langevin MD will be ever more significant as the number of particles increases, both at and away 
from the critical point. 

7 Discussion 

We have described a Fourier-based Langevin scheme, capable of accelerating the dynamics of par- 
ticulate systems with a Lagrangian data representation. We have demonstrated that there is great 
potential in speeding up the dynamics of long-wavelength modes. Two issues related to the ac- 
celerated dynamics will be addressed in future research. (1) Dynamical faithfulness: how much of 
the actual dynamics is unchanged by taking different time steps for different wavelength modes? 
(2) Does this method offer even more of a gain when it is applied to molecular systems with 
nonconserved order parameters (for example, dipolar systems or systems with structural phase 
transitions)? We believe that our preliminary computational investigations Fourier methods have 
shown great promise, and we intend to explore these questions in detail in future work. 
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